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Abstract — We consider the problem of designing (or augment- 
ing) an electric power system at a minimum cost such that 
it satisfies the N-k-E survivability criterion. This survivability 
criterion is a generalization of the well-known N-k criterion, and 
it requires that at least (1 —£j) fraction of the total demand 
to be met after failures of up to j components, for j = 1 , ■ • • , k. 
The network design problem adds another level of complexity 
to the notoriously hard contingency analysis problem, since 
the contingency analysis is only one of the requirements for 
the design optimization problem. We present a mixed-integer 
programming formulation of this problem that takes into account 
both transmission and generation expansion. We propose an al- 
gorithm that can avoid combinatorial explosion in the number of 
contingencies, by seeking vulnerabilities in intermediary solutions 
and constraining the design space accordingly. Our approach 
is built on our ability to identify such system vulnerabilities 
quickly. Our empirical studies on modified instances from the 
IEEE 30-bus and IEEE 57-bus systems show the effectiveness 
of our methods. We were able to solve the transmission and 
generation expansion problems for k = 4 under 2 minutes, while 
other approaches failed to provide a solution at the end of 2 
hours. 

Index Terms — Long-term grid planning, contingency require- 
ments, decomposition, separation oracle, implicit optimization. 



I. Introduction 

CONTINGENCY analysis of power systems is a problem 
of increasing importance due to ever increasing demand 
for electricity. This increase in elecnicity demand outpaces the 
increase in system capacity growth, which leaves less room 
for redundancy in the system. As a consequence, equipment 
failures are now more likely to lead to blackouts. At the 
same time, society's increased reliance on electricity breeds 
far and wide ramifications for a blackout. While minimizing 
the likelihood of any blackout is a laudable goal, this will 
inevitably come at an economic cost. A partial remedy for 
increasing costs is to take into account both the likelihood 
of an event and severity of its consequences in the planning 
phase. This concept is adopted in Transmission Planning 
Standard (TPL-001-1, [1]), defined by the North American 
Electric Reliability Corporation (NERC). According to this 
standard, if only a single element is lost (N — 1 contingency), 
the system must be stable, without any loss-of-load. In the case 
of k simultaneous failures {N — k contingency), the system still 
has to restore stability, but a limited loss-of-load is allowed. 
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While contingency analysis, as a concept, is simple, al- 
gorithmically it is far from trivial. Despite its difficulty, the 
critical importance of the problem has drawn a lot of interest, 
and significant progress has been made over the last decade. In 
particular, optimization methods have been proposed to replace 
enumerative approaches that rely on verifying feasibility on 
each state of a given list of contingencies. For example, 
optimization methods are used to find small groups of lines 
whose failure can cause a severe blackout or a large loss- 
of-load [2], [3], [4], [5], [6]. N-k contingencies were also 
studied in optimal power flow models [7], [8], [9] and unit 
commitment problems [10] using a single bus model. The 
methods used in [6], [9], [10] are all based on a bilevel 
programming approach, which is the main method used for 
network inhibition/interdiction problems. Probability analysis 
[11], limitations on generation and line capacity [12], swarm 
optimization [13], and methods based on topological charac- 
teristics of power grids [14] have been used for vulnerability 
analysis. Other approaches for contingency analysis are re- 
viewed in [15], [16]. 

The great strides in contingency analysis over the last 
decade have established the basis for higher objectives. In this 
paper, we consider the transmission and generation expansion 
planning (TGEP) problem with contingency constraints. Our 
goal is to design (or augment) a power system at a minimum 
cost to satisfy contingency constraints. We formalize the con- 
tingency requirements with the N-k-e criterion, where e is a 
parameter vector that specifies allowable loss-of-load, for each 
contingency size, as a fraction of total system demand. More 
specifically, this criterion requires that for any contingency of 
size J = 0, lv ,k, at least (1— £,■) fraction of the total demand 
must be satisfied. Following NERC's TPL-001-1 standards, for 
the no-contingency state and contingencies of size one (j = 1), 
no loss-of-load is allowed (i.e. £o = £\ = 0); for multiple failure 
contingencies (j > 2), a small fraction of total load demand 
can be shed (i.e. < £2 < • • • < £k < 1)- 

To understand the complexity of the TGEP, one should ob- 
serve that contingency analysis, a difficult combinatorial opti- 
mization problem by itself, is only one of the prerequisite steps 
in solving TGEP; we must also address network design issues. 
Our approach is built on our ability to solve the contingency 
analysis problem efficiently so that it can be embedded in a 
broader framework. We begin by formulating a mixed-integer 
linear program (MILP) to model TGEP that explicitly includes 
multiple states representing each of the possible contingency 
states and, for each of these contingency states, contains the 
corresponding generation and flow variables to ensure that 
at least (1 — £;) fraction of the demand can be met. Such 
an enumerative approach, however, becomes computationally 
intractable even for moderate values of system size N and 
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maximum contingency-size budget k, since the number of 
states grows with (v). Furthermore, it should be noted that 
contingency analysis for TGEP requires looking at higher 
values of k, compared with contingency analysis for shorter 
term problems. The reason is that in long term problems, 
planned outages, such as system components going under 
maintenance, contribute to k, whereas in shorter-term problems 
they do not. For instance, for day-ahead planning, we will 
know the unavailable system components in advance; remove 
them from the system; and then we will do contingency analy- 
sis for unexpected failures in the remaining system. However, 
for long term problems such as TGEP, we want to design a 
system such that a feasible direct current optimal power flow 
(DCOPF) exists for any combination of planned and unplanned 
outages. This forces us to look at higher values of k, which 
makes enumerative approaches prohibitively expensive, even 
when state of the art high performance computing platforms 
are available. 

To overcome this challenge, we propose two cutting plane 
algorithms, one based on a direct application of Benders 
decomposition method to check the load satisfaction of each 
contingency state explicitly and another based on an online 
contingency state generation (OCS) algorithm, which solves 
bilevel separation problems to determine the worst-case loss- 
of-load for each contingency size j = 0, 1, • • • ,k. The OCS 
algorithm implicitly identifies worst-case contingencies with- 
out explicitly evaluating all contingency states, and generates 
additional constraints, corresponding to Benders feasibility 
constraints, to exclude solutions that are infeasible under the 
identified worst-case contingencies. This approach avoids the 
prohibitive cost of generating a constraint for every combina- 
tion of failures, by only identifying relevant combinations as 
needed. As the computational experiments show, this reduction 
to only relevant contingencies makes a tremendous impact on 
the scalability of the algorithm. 

We applied these approaches to the IEEE 30-bus and IEEE 
57-bus systems. Computational results show the scalability of 
the OCS algorithm. We were able to solve the TGEP problem 
for up to k = 4 in under 2 minutes, while explicit enumeration 
approaches, e.g. extensive form and Benders decomposition, 
failed to provide any solutions at the end of 2 hours. We 
observed the key to our success was our ability to avoid 
looking at all individual vulnerabilities, and only a small 
number of iterations (vulnerability searches) are enough to 
find a provably-optimal solution. 

Transmission and generation expansion problems have been 
studied for a long time. A survey of earlier work in this 
area can be found in [17]. Recently, especially after the 
2003 Northeast American blackout, security, stability and 
reliability issues have been becoming another major concerns 
in modern power systems. Our earlier results were presented 
in [18], [19]. [20] and [21] addressed the defense protection 
of power system, and they analyzed the interaction between 
a power system defender and a terrorist who seeks to disrupt 
system operations. [22], [23], and [24] studied the contingency 
criteria by stochastic programming and integer programming 
approaches. [25] and [21] proposed a multilevel mixed integer 
programming models for transmission expansion. Most of 



the work in this area however, is restricted to the N — I 
contingency criteria, which cannot reflect current situation for 
contingencies with more than one failure. Additionally, the 
proposed models cannot be directly extended for this new 
situation, and they only consider the transmission expansion, 
and contingencies are restricted to transmission elements. 

The generation expansion problem has recently been studied 
in [26]. Transmission expansion and generation expansion 
planning problems have also been studied in a unified model 
[27]. Integration of renewable energy resources was taken into 
account generation and transmission expansion planning [28], 
[27], [26]. However, none of these studies consider the contin- 
gency criteria in case of failures of both generating unit and 
transmission elements. 

The rest of this paper is organized as follows. In Section 
II, the TGEP problem considering the full set of contingency 
states is formulated as a mixed integer nonlinear program 
(MINLP); Section III presents two methods to solve this large- 
scale MINLP; in Section IV, results of numerical experiments 
performed on two IEEE test systems are presented; Section V 
concludes the paper. 

II. Models 

A. Nomenclature 
Sets and indices 

V Set of buses (indexed by i,f). 

Sj Set of all contingency states with exactly j failures 

(indexed by s). 
S Set of all contingency states with k or less failures, 

S = S 1 l)S 2 ---US k . 
\s\ Number of failed element(s) in contingency s. 
G Set of generating units (indexed by g). 
Gj Set of generating units at bus i. 
E Set of transmission elements (indexed by e). 
E.i Set of transmission elements oriented into bus i. 
Ei. Set of transmission elements oriented out of bus ;. 
i e ,j e Tail/head (bus no.) of transmission element e = (i e ,j e ). 

Parameters 

C e Investment cost of transmission element e. 

C g Investment cost of generating unit g. 

Cg Marginal production cost of generating unit g. 

P g Maximum capacity of generating unit g. 

B e Electrical susceptance of transmission element e. 

F e Capacity of transmission element e. 

Di Electricity load demand at bus i. 

D Total electricity load demand across all bus ; G V, 

D = ZievDi. 
M e Big M value for transmission element e. 
a Weighting factor to make investment cost and 

operating cost comparable. 
N Total number of transmission elements and 

generating units, N = \G\ + \E\. 
k Maximum contingency size under consideration. 

In any contingency state the number of failures, i.e. the 

number of system elements in a contingency, is 
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between and k. 
Ej Fraction of load demand that can be shed given size j 

contingencies, for all j = 0, 1, • • • ,k. 
d s g Binary parameter that is 1 if generating unit g is 

part of the contingency state s and otherwise. 
d s e Binary parameter that is 1 if transmission element e 

is part of the contingency state s and otherwise, 
d Vector that concatenates d s g and d s e variables. 

Decision variables 

x g Binary variable that is 1 if generating unit g is added 

and otherwise. 
x e Binary variable that is 1 if transmission element e is 

added and otherwise, 
x Vector that concatenates x g and x e variables. 
f e Power flow on transmission element e. 
p g Power output of generating unit g. 
qi Loss-of-load at bus /'. 
9i Phase angle of bus i. 

d g Binary variable that is 1 if generating unit g is 

selected to be in the contingency and otherwise. 
d e Binary variable that is 1 if transmission element e is 

selected to be in the contingency and otherwise. 
f s e Power flow for transmission element e for 

contingency state s. 
p g Power output of generating unit g for contingency 

state s. 

q\ Loss of load at bus i for contingency state s. 
6j Phase angle of bus i for contingency state s. 

For a contingency state s G S, d* = 1 and d s e = 1 denote that 
generating unit g and transmission element e fail in this state, 
respectively. Conversely, d g — and d s e = denote that both 
these two elements are available. Thus, for 5 = denoting the 
no-contingency state, we have "5 = and d® = for all g € G 
and e G E. 



B. Transmission and Generation Expansion Model 

In this section, we extend the standard TGEP problem to 
include contingency constraints. For brevity of presentation, 
we treat all power system elements, including existing ones, as 
candidates for addition to the system. For an existing element, 
the investment cost (C e ,C g ) is set to and the corresponding 
investment decisions (x e ,x g ) fixed at 1. 

Once the network design decisions are made, each element 
selected for addition becomes available in all contingency 
states s G Sj and j = 0, 1, • • ■ ,k, unless it is part of a given 
contingency. Consistent with the NERC reliability standards, 
in the no-contingency state {N-0) and the single contingency 
state (AM), no loss-of-load is allowed, thus £o = £i = 0. For 
larger contingency states, e.g. j =2, ■■■ ,k, the total loss-of- 
load is limited by the threshold £,-. We assume that £j < £; + \ 
for all j = 0,1,--- ,k- 1. The MINLP model for TGEP is 
formulated as follows, 



min V C e x e + V CgXg la[ft! (la) 
x,f,p,q,8 £ E g % g ^ G 

SX - E Pl + £ fe ~ I fe + <fi = D U Vi G V, (lb) 

gGG, e&Ej eeEi, 

seS 

B e (6l-6 s k )x e {\-d s e )-f e =Q, "ie EE, (lc) 

seS 

-F e x e {l-d s e )<f e <F e x e {\-d s e ), VeG/s, (Id) 

seS 

Q<Pg<P g x g (\~d s g ), VgeG,seS (le) 

£ < £|. S |D, VseS (If) 
iev 

0<q S i<D h VieV,\/seS (lg) 

*,e{o,i}, VgGG (ih) 

x e e{0,l}, VeeE (li) 

In all subsequent formulations, unless otherwise specified, the 
indices i,g,e,s and k are elements of their corresponding sets, 
i.e., ;' G V,g G G,e G E, and s G S. 

The objective (la) is to minimize the total transmission and 
generation investment cost plus the weighted operating cost 
in the no-contingency state (s = 0). Note that our formula- 
tion does not take into account the operational costs during 
a contingency state. There are two reasons for this. First, 
contingencies have low likelihood and thus the operational 
costs during such events are negligible. The real financial 
burden of a contingency shows itself when the contingency 
leads to a blackout, which brings us to the second reason: 
the primary goal during a contingency is to keep the system 
intact as opposed to minimizing operational costs, since the 
cost of system failure is likely to be significantly higher than 
any operational cost. Constraints (lb) are conservation of 
flow requirements for each bus and contingency pair. For any 
transmission element that is operational, Kirchhoff's voltage 
law must be enforced by (lc). Power flow on transmission 
element e is governed by thermal capacity constraints (Id). For 
each contingency state, the power output of a generating unit 
must satisfy the upper bound given by (le). Constraints (If) 
dictate that the total loss-of-load in the system cannot exceed 
£| s | fraction of the total system load, i.e., at least (1 — £u)D 
of demand must be satisfied. Constraints (lg) restrict the loss- 
of-load at each bus to be at most the total demand at the bus. 

Observe that constraints (lb)-(lg) are specific to a particular 
contingency state; that is, for a given contingency state s, 
the transmission and generation element(s) in the contingency 
have zero capacity. In the no-contingency state s = 0, all in- 
vested transmission elements and generating units are available 
for the DCOPF problem. 

III. Solution Approaches 
Replacing constraints (lc) by 

B e (d s ie -d s je )-f e +M e {\-x e + d s e )>Q, Ve,s, (2) 
B e {ef e -e s ^-f e -M e {l-x e + d s e )<Q, Ve,s. (3) 
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where M e is a sufficiently large constant, transforms formu- 
lation (1) to a mixed integer linear program (MILP), which 
we refer to as the extensive form (EF). EF will typically 
have an extremely large number of variables and constraints 
because it grows with the number of contingency states, which 
increases exponentially with N and k. For large power systems 
and/or a contingency k greater than one, EF rapidly becomes 
computationally intractable for increasing system size, N, and 
increasing contingency size, k. In the following sections, we 
modify this formulation and present cutting plane algorithms 
for solving the reformulated problem. 

A. Benders Decomposition 

We begin by presenting an alternative formulation with only 
|G| + \E\ binary variables but possibly an extremely large 
number of constraints. We use linear programming duality to 
generate valid inequalities for the projection of the natural 
formulation onto the space of the x variables. In essence, we 
use a variant of Benders Decomposition, in which we only 
generate valid inequalities corresponding to "feasibility" cuts. 

For a given contingency state, s G S, capacity expansion vec- 
tor x and contingency state vector d\ we solve the following 
linear program, denoted as the primal subproblem PSP(x,d i ), 
to determine a DCOPF that minimizes total system loss-of- 
load. 



min £ q\ 

f,p,q,9 



(4a) 



s.t. (of) £ p g + £ ft - £ ft + q\ = Dj, V/ (4b) 

geG/ eeE; eeE L 

(#) -B e (8? e -e s je )+ft<M e (l-x e +d s e ), Ve (4c) 

(#) B e (8? e -8l)-ft<M e (l-x e + d s e ), Ve (4d) 

(5/) ft<F e x e (l-d s e ),Ve (4e) 

(rtf) -ft<F e x e (l-d s e ), Ve (4f) 

(Q) 0<p s g <P g x g (l-d s g ),\fg (4g) 

W ) < q\ < D t , Vi (4h) 

In this formulation, the objective, (4a), is to minimize total 
loss-of-load by adjusting the flow, phase angles and power 
generation, given the prescribed capacity expansion decision x 
and contingency vector d s , corresponding to scenario s. Letting 
z(x,d s ) be the optimal objective value of (4), if z(x,d s ) > £\ S \D, 
there does not exist a feasible DC power flow satisfying at 
least (1 — £|i|) of total demand. Alternatively, if z(x,d s ) < £\ S \D, 
there exists a power flow that can satisfy at least (1 — £u) of 
total demand. 

Variables in parentheses on the left-hand-side of the con- 
straints in (4) denote the corresponding dual variables. In turn, 
we can formulate the dual of this problem, DSP(x,cF) as 
follows, 

max £A(of + 1^(1 -Xe + d s M + fc) 

a, ji.p.S, 7j,f,A ieV eeE 

+ £ F e x e (\ - d s e w e + n s e ) + £ p gXg (i - d s g )Q, 

eeE g£G 



subject to constraints corresponding to primal variables 
f,p,q, 9. Since PSP(x,d s ) has a finite optimal solution value 
(in the worst case, all load will be shed), DSP(x,d v ) also 
has a finite optimal solution value and by strong duality, 
the optimal solutions coincide. Since DSP(x,d s ) has a finite 
optimal solution value it also has an optimal extreme point. 
Thus we can reformulate PSP(x,d v ) as follows, 

max £/),(«/ + A/) + £M e (l -x e + d s M + ft) 

leL ieV e£E 

+ £ F e x e (l-dM + ni) + £ P g x g ( 1 - ^) $ , (5) 

eeE g£G 

where L" is the set of extreme points corresponding to the 
polyhedron characterized by dual constraints based on (4) for 
primal variables f,p,q,0. 

The constraint z(x,d s ) < £\ S \D should be satisfied for all 
s G S. Thus contingency feasibility conditions can be defined 
as follows, 

£ Di (af + A/)+ £m,(1 -x e + d°M + p?) 

i€V eeE 

+ £ F eXe {\ - d s M + n l) + £ p gXg {\ - (6) 

eeE geG 

<e ]s] D, \/eeL s ,seS 

Since the objective is to minimize total investment cost and 
operating cost in the no-contingency state (s = 0), we enforce 
all constraints for state explicitly. For all other contingency 
states s <E S\ {0}, constraint set (6) ensures that at least 
fraction of the total demand is satisfied. The reformulation of 
(1) is given as: 



min £ C e x e + £ CgXg + a £ q>° (7a) 

*,f,P,q,8 e tk geG gTc 

s.t. £D,-(«/ + A/) + £M e (l + 

iev eeE 

+ £ F e x e {\ aZ){8l + nl) + £ P g x g (l ~ d s X 

eeE geG 

<£| S |D, VeeL s ,seS\{0} (7b) 
L^+L/^L/^ V/6V (7c) 

geGj eeEj eeEj, 

B e (ef e - 8l) -tf+M e {\ -x e ) > 0, Ve g E (7d) 
B e {el - 0l) -f)-M e {\ -x e ) < 0, Ve G E (7e) 



-F e <f e <F e , VeeE 
0<p° s <Pg X g, VgeG 
x g G{0,l}, VgeG 
x e G {0,1}, Ve EE 



(7f) 
(7g) 
(7h) 
(7i) 



The number of constraints in formulation (1) grows expo- 
nentially with the problem size, so we solve it via Benders 
Decomposition (BD). At a typical iteration of BD, we consider 
the restricted master problem (RMP) (7), which has the same 
objective as (1) but involves only a small subset of the 
constraints in (7b). We briefly outline BD below. For a detailed 
treatment of BD please refer to [29]. 
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Let t be the iteration number and let the initial RMP 
be problem (7) without any (7b) constraints. Let x' be a 
concatenation of the expansion variables in the f th iteration. 



Algorithm 1 Benders Decomposition (BD) 

l: f<-0 

2: solve RMP 

3: if RMP is infeasible 
4: EXIT, TGEP is infeasible 

5: else 

6: Let x f be the optimal solution to RMP 

7: for s £ S, solve DSP^er 5 ), let z' be the objective value 

8: if z' > S\ S \D 

9: Add feasibility cut (6) to RMP 

10: end if 
11: end for 

12: if no feasibility cut(s) added in Step 9 
13: x' is optimal, EXIT 
14: else 

15: t <-t + l, go to Step 2 
16: end if 

17: end if 



By using a Benders reformulation, we are able to decom- 
pose the extremely large formulation (1) into a master problem 
and multiple subproblems (one for each contingency state). 
In theory, this enables us to solve larger instances, which 
would not be possible by a direct solution of EE However, 
the extremely large number of contingency states makes direct 
application of Benders ineffective for large power systems 
and/or a non-trivial contingency budget (i.e., k > 1). In the 
next section, we develop a custom cutting plane algorithm 
that evaluates all possible contingency states implicitly using 
a bilevel separation oracle. 

B. Online Contingency Scenario Generation 

Sizes of most power systems in operation (typically thou- 
sands of generating units and transmission elements) may 
preclude direct solution of (1). Even using a decomposi- 
tion algorithm (e.g. BD) may not be feasible because each 
contingency state must be considered explicitly. Our goal is 
to instead use a separation oracle that implicitly evaluates 
all contingency states and either identifies a violated one (a 
contingency with j failures (for all j = 1 , • • • , k) that cannot 
be survived by the current power system design) or provides 
a certificate that no such contingency state exists. If such 
a contingency exists, we use this contingency to generate a 
violated Benders feasibility cut, as described in the previous 
section, for the RMP. If no such contingency exists, then the 
current capacity expansion x is optimal and we terminate the 
algorithm. 

1) Power Sy stern Inhibition Problem (PSIP): Given a ca- 
pacity expansion decision x, the Power System Inhibition 
Problem (PSIP) can be used to determine the worst-case 
loss-of-load under any contingency with j failures, for all 
j = !,■■■ ,k. In this bilevel program, the upper level decisions 
d correspond to binary contingency selection decisions and the 



lower level decisions (f,p,q,0) correspond to recourse power 
flow, generation scheduling, and load shedding decisions rel- 
ative to the given contingency, prescribed by d. 

Note that in the prior model d J was an input parameter, 
whereas in this formulation, we are now selecting the elements 
of the contingency, with d becoming a vector of decision vari- 
ables. For clarity of exposition, the superscript s corresponding 
to variables f , p, q and has been removed, as the contingency 
state is not pre-specified, but rather part of the decision making 
process within the PSIP(x, j). PSIP(x,y) is given as follows: 

max min } a, (8a) 

d f,p,q,6 

s.t. E d e + E dg = j, (8b) 

e£E geG 

(«<) E Pg + E fe - E fe + 91 = A. V/ (8c) 

geGj eeEj ee£, : 

(fie) - B e ( 9 ie - d je )+f e <M e (l-X e + d e ),Ve (8d) 
(fie) B e (9 ie - 9 je ) -f e < M e {\ -x e + d e ), Ve (8e) 



(Se) f e <F e x e (l-d e ),Ve (8f) 

(rje) -f e <F e x e (l-d e ),Ve (8g) 

(&) 0<p g <P g X g (l-d g ), Vg (8h) 

(Xi) 0<qi<Di, Vz (8i) 



The objective (8a) is to maximize the minimum loss-of-load. 
For a given contingency state defined by d, the objective of the 
power system operator (the inner minimization problem) is to 
determine the DCOPF such that the loss-of-load is minimized. 
Constraint (8b) is a budget constraint limiting the number 
of power system elements that can be in the contingency. 
Constraints (8c) are standard flow conservation constraints. 
Constraints (8d) and (8e) together enforce Kirchhoff's voltage 
law, for active transmission elements. Constraints (8f) and 
(8g) are constraints associated with the capacity of each 
transmission element. Constraints (8h) limit the maximum 
capacity of each generating unit. If a generating unit g is 
NOT part of the contingency (i.e., d g = 0), then the maximum 
capacity of the generating unit is enforced, if the unit was 
added (x g = 1). Otherwise, the power output of the generating 
unit must be zero. 

The upper-level decisions of this bilevel program are to se- 
lect a contingency, using the binary variables d, that maximizes 
the subsequent loss-of-load in the lower-level problem. 

2) A MILP Reformulation of PSIP: Bilevel programs like 
(8) cannot be solved directly. One approach is to reformulate 
the bilevel program by dualizing the inner minimization prob- 
lem. For fixed values of d, the inner minimization problem is 
a linear program that is always feasible. By strong duality, and 
combining the upper level of (8) we can reformulate the bilevel 
program as a single level bilinear program, where the objective 
function of the dualized problem contain terms associated 
with the product of the upper-level contingency selection 
variables d and the lower-level dual variables (j3,/3,(5,77,£) 
associated with flow balance, transmission flow, transmission 
capacity, and generation capacity constraints. However, with 
additional variables and constraints, these bilinear terms can 
be linearized. 
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Each bilinear term can be linearized using the following 
strategy. Let u < and v < be continuous variables and b € 
{0, 1}. Then the bilinear term, bu, can be linearized as follows. 
Letting v = bu, we introduce the following three constraints to 
linearize the bilinear term bu. 



v > u-U(l-b) 

v > -Ub 

v < u + U(l-b) 



(9a) 
(9b) 
(9c) 



Here, parameter U represents a valid upper bound for con- 
tinuous variable u and satisfies U > \u\. Assessing these three 
constraints for both binary values of b show that they indeed 
provide a valid linearization. If b = 0, then constraints (9b) 
and v < together imply that v = 0. With v = 0, constraints 
(9a) and (9c) together imply that — U <u<U, which are never 
binding. If b = 1, then constraints (9a) and (9c) together imply 
u = v and constraints (9b) and v < implies —U < v < 0, which 
is valid. 

We follow a similar strategy to linearize all five bilinear 
terms (j5d e ,j5d ei 8d e ,rjd e ,^dg). Define continuous variables 
(r\r 2 ,r 3 ,r 4 ,r 5 ) and let r\ = $ e d e , r 2 e = $ e d e , r\ = 8 e d e , r 4 e = 
rid e , and r^ = t, g d g . Following the same linearization strategy 
introduced in (9), we now state the full mixed integer linear 
PSIP formulation for completeness, which we call the Mixed- 
Integer Power System Inhibition Problem, M-PSIP(x,y): 

max !>,(«,■ +Xi) + £ M e {\ -X e ){Pe + Pe) 

d.a.^p.S.n^.r i e v eeE 

+ L (M e {r\ + n) + F e x e (S e + vie) - F e x e (r 3 e + r 4 e ) 

eeE V 

gee 

s.t. £<4+£4=./> (10b) 

eeE geG 

ccj e - a ie +p e -p e + 8 e -ri e = 0, Me (10c) 

Oi g + Cg<0, Wg (lOd) 

Ck + Xi < 1, V/ (lOe) 

£ B e (pe-Pe)+ £ B e (P e -P e )=0,Vi (lOf) 

Pe < 0, j3 e < 0, 8 e < 0, r, e < 0, Ve (lOg) 

r\ <0, r] <0, rl < 0, r 4 < 0, Ve (lOh) 

C A -<0, r e 5 <0, Vg (101) 

h < 0, V/ (lOj) 

Next, we outline an algorithm for optimally solving problem 
(1) that combines a Benders decomposition with the aid of an 
oracle given by (10), which acts as a separation problem. A 
given capacity expansion x is optimal if the oracle cannot find 
a contingency of size j, for any j = 1 , • • • , k, that results in a 
loss-of-load above the allowable threshold EjD. 

For each contingency budget j, we can check for y'-element 
contingencies by solving M-PSIP using a failure budget of 
j (i.e. the right-hand side of inequality (10b) is set to j). 
Whenever the oracle determines that the capacity expansion 
decision x is not N-k-e compliant, it returns a contingency 



d that results in a loss-of-load, above the allowable threshold 
EjD for y'-element failures. 

Let r be the iteration number and let the initial RMP 
be problem (7) without any (7b) constraints. Let x r be a 
concatenation of the expansion variables (x r e ,x r g ). 

Algorithm 2 Online Contingency Screening (OCS) 
t^0 

Solve RMP for iteration t 
if RMP is infeasible 

EXIT, TGEP is infeasible 
else 

Let x' be the optimal solution to RMP 
for j = !,-■■ ,k 

Solve M-PSIP(x f ,;'), let z'j be the objective value 
and d'j be the contingency selection decision 
if z?j > EjD then 

Solve DSP(x',<r}), add feasibility cut (6) to RMP 
t 4— t -h 1, go to step 2 
end if 
end for 
end if 

x ? is optimal, EXIT 



At each iteration, either a contingency that results in loss- 
of-load above the allowable threshold is identified and a 
corresponding feasibility cut is generated and added to RMP, 
or no contingencies are found, which means that the current 
solution is optimal and the algorithm can terminate. 

IV. Numerical Experiments 

We implemented the proposed models and algorithms in 
C++ and CPLEX 12.1 via ILOG Concert Technology 2.9. All 
experiments were run on a machine with four quad-core 2.93G 
Xeon with 96G of memory. For the following computational 
experiments, a single CPU and up to 8GB of RAM were 
allocated. The optimality gap was set to be 0.1% for CPLEX. 

We have tested our models and algorithms on the IEEE 30- 
bus and IEEE 57-bus systems [30]. For each power system, 
we consider five different contingency budgets k = 0,1,2,3, 
and 4. Altogether, we consider 10 instances. 

Table I compares the run times for the three different 
approaches. For each of the 10 instances, m provides the 
number of distinct contingencies. Initially for each test sys- 
tem, we replicate a subset of existing generating units and 
transmission lines to create a set of candidate elements. With 
these candidate elements as a starting point, we iteratively 
solve the PSIP problem for k and e values presented in 
Table I using OCS. In all 10 instances, the values of £, are 
invariant. That is, for instance, &i = 0.05 when k was chosen 
to be 2,3, or 4. For these, we identify vulnerabilities in the 
power system and introduce additional candidate generation 
and transmission elements. We follow this method to create 
the augmented the IEEE 30-bus and IEEE 57-bus test systems 
for the computational experiments presented subsequently. 
We want to note that our contributions in this paper are in 
algorithmic fundamentals, and thus are not very sensitive to 
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the particular problem instances. Hence, the particular list of 
candidates lines for a problem instance will only have a minor 
effect on the performances of our algorithms. 



TABLE I 

Run times for different solution approaches 



the current network design; and the dual subproblems (DSP), 
which generates the feasibility cuts. 



Test Systems 




k 


e 


EF 


BD 


ocs 


IEEE 30-bus 






















152 


1 





76 


2 


4 




> UK 


2 


0.05 


X 


87 


16 




> 500AT 


3 


0.10 


X 


3,126 


128 




>21M 


4 


0.20 


X 


X 


141 


IEEE 57-bus 






















110 


1 





27 


164 


120 




>5K 


2 


0.05 


X 


870 


36 




>200K 


3 


0.10 


X 


X 


51 




>5M 


4 


0.20 


X 


X 


65 



Solution time (sees) 



Table I provides the run time (in CPU seconds) for each 
instance under the three different approaches. In this table "x" 
means the algorithm failed to complete at the end of 2 hours. 
Note that the first approach, the extensive form (EF), can only 
solve the smallest of instances. This is because of the sheer 
size of the problem, in which, for each contingency, a full 
DCOPF problem must be embedded in the formulation. As 
the number of contingencies grows, this formulation quickly 
becomes intractable. Note that we are working with small data 
sets here and target systems will be in the order of thousands 
of elements, which will make EF intractable even sooner. 

The second approach, BD, bypasses this problem via a 
Benders decomposition, with corresponding delayed cut gener- 
ation. However, this still suffers from the combinatorial growth 
in the number of contingency states - for each contingency, 
a subproblem (DSP) must be solved to check for violated 
feasibility cuts to add to the RMP We see that larger problem 
instances can be solved, relative to EF, but the BD approach 
nonetheless cannot solve the largest problem instances. 

With the OCS approach, we see that all instances of the 
problem can be solved, in all cases in under three minutes 
and frequently in only a few seconds. This is a result of the 
combination of the strength of the Benders cuts, enabling the 
problem to be solved in a very limited number of iterations, 
and also the fact that we are able to implicitly evaluate the 
contingencies in order to identify a violated contingency and 
then quickly find its corresponding feasibility cut by solving 
a single linear program (DSP). 

Table II provides us with further evidence for the scalability 
of OCS. For each instance, we see the total number of possible 
contingency states m and then the number of contingency 
states for which corresponding feasibility cuts were actually 
generated, denoted by 'cont.' in the table. Clearly, it is a 
very tiny fraction of the possible number of contingencies, 
which is critical to the tractability of the approach. The 
remaining columns of this table breakdown the total run time 
by time spent on the three components of the algorithm - the 
RMP, which identifies a candidate network design; the mixed- 
integer linear power system inhibition problem (M-PSIP), 
which identifies a contingency that cannot be overcome by 



TABLE II 
OCS RUNTIME BREAKDOWN 



Test Systems 


m 


k 


e 


RMP 


M-PSIP 


DSP 


cont. 


IEEE 30-bus 

























152 


1 





1 


3 





2 




> UK 


2 


0.05 


1 


15 





5 




>50QK 


3 


0.10 


1 


127 





8 




> 21M 


4 


0.20 


1 


140 





10 


IEEE 57-bus 

























110 


1 





83 


37 





3 




>5a: 


2 


0.05 


15 


21 





4 




>200K 


3 


0.10 


14 


37 





7 




>5M 


4 


0.20 


14 


51 
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V. Conclusion 

We studied the transmission and generation expansion prob- 
lem with contingency constraints. More specifically, we inves- 
tigated the problem of improving an electric power system at 
a minimum cost by adding generators and transmission lines, 
such that it satisfies the N-k-e survivability criterion. This 
survivability criterion is a generalization of the well-known N- 
k criterion, and it requires that at least (1 — £j) fraction of the 
total demand is met even after failures of any j system com- 
ponents, for all j = !,-■■ ,k. This design problem adds another 
level of complexity to the contingency analysis problem, since 
the contingency analysis is only one of the constraints in the 
design optimization problem. We proposed two algorithms: 
one is based on the Benders decomposition approach, and 
the other is based on online contingency state generation. 
The latter approach avoids the combinatorial explosion by 
seeking vulnerabilities in the current solution, and generating 
constraints to exclude such infeasible solutions. We tested our 
proposed approaches on the IEEE 30-bus and the IEEE 57- 
bus systems. Computational results show the proposed online 
contingency generation algorithm, which uses a bilevel sepa- 
ration technique to implicitly consider all exponential number 
of contingencies, significantly outperforms a standard Benders 
decomposition. We were able to solve all instances in our 
experiment in under three minutes, while the extensive form 
and the Benders Decomposition algorithm failed to complete 
at the end of 2 hours. 

We believe that this paper will provide the fundamentals 
for many other studies in contingency-aware transmission and 
generation expansion. As an example, we want to apply for 
methods to full-scale systems. While our results are very 
promising in terms of scalability, full-scale problems will 
surely pose some computational challenges and will require 
adopting high-performance computing resources. Also our 
current model assumes all failures happen simultaneously. In 
order to reflect practical operation situations, where failures 
may happen consecutively, new models that consider timing 
between system element failures are needed. Additionally, 
unit commitment and de-commitment is not considered in 
our current model. We plan to extend our models for these 
cases. Finally, we worked with a deterministic model, and 
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it is essential, to take stochasticity into account for planning 
problems. We believe our current frame work can be naturally 
extended for stochastic problems. 
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